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Abstract 



This paper tests a dielectric model for variation of hydration free energy with 

in 

Q"^ ' geometry of complex solutes in water. It works out some basic aspects of the 



theory of boundary integral methods for these problems. One aspect of the 
algorithmic discussion lays the basis for multigrid methods of solution, meth- 
■ ods that are likely to be necessary for similarly accurate numerical solution of 

these models for much larger solutes. Other aspects of the algorithmic work 
^ • show how macroscopic surfaces such as solution interfaces and membranes 

may be incorporated and also show how these methods can be transferred 
directly to periodic boundary conditions. This dielectric model is found to 
give interesting and helpful results for the variation in solvation free energy 
with solute geometry. However, it typically significantly over-stabilizes clas- 
sic attractive ion-pairing configurations. On the basis of the examples and 
algorithmic considerations, we make some observations about extension of 
this continuum model incrementally to reintroduce molecular detail of the 
solvation structure. 
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1. Introduction 

This paper tests a dielectric model for variation of hydration free energy with geom- 
etry of complex solutes in water. It works out some basic aspects of the theory of boundary 
integral methods for these problems and describes, demonstrates, and tests such algorithms 
for solution of that macroscopic Poisson equation. It then makes some observations about 
extension of this continuum model incrementally to reintroduce molecular detail of the sol- 
vation structure. 

Numerical solution of these models has been discussed many times before 0-25|. Di- 



electric continuum models are widely considered simplistic descriptions of the solvation of 
charged or polar solutes in solutions. Thus, it might not be obvious whether accurate so- 
lutions to the governing Poisson equation are important. However, when solved to thermal 
accuracy these models can give helpful and interesting descriptions of the variation of solva- 
tion free energy with geometry p6ip7| . These dielectric models are physical. Furthermore, 
their connection with molecular theory is known and straightforward p6 H28|j3TH 34| . For 
these reasons the dielectric models deserve to be exploited more thoroughly than they have 
been. These interesting properties of the physical model require methods that yield solvation 
free energies accurate on the thermal energy scale (< ksT ^ 0.03e\^ under the conditions of 
widest interest). Experience has shown that this level of accuracy can be difficult to achieve 



and requires careful numerical methods pI| , |23| -[25[] 



The Poisson equation is a linear equation and the chief difficulty is associated with the 
treatment of interface boundary conditions on a molecular surface that can be irregular. 
The possibility of disconnection of the molecular surface with molecular rearrangement is of 
utmost physical importance. Viewed from a physical level, those disconnections can produce 
microscopic solvent pockets that the model treats as bulk continuum solvent. It might be 
questioned whether the model is significantly less accurate in those circumstances. Viewed 
on a technical level, those disconnections greatly complicate the description of the molecular- 
solvent interface. Because spatial resolution in the description of that interface is important 
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to achieving accuracy in the solutions, it might be further questioned whether it is feasible or 
desirable to obtain solutions of the models that are sufficiently accurate to examine the effects 
of those disconnections. But boundary element approaches are expected to be particularly 
appropriate for these purposes. This has been clearly expressed by Yoon and Lenhoff ||14|| . 
The development here follows that idea but accomplishes the required boundary (or surface) 
integrations utilizing two simple procedures: (a) plaques are defined on the basis of quasi- 
random number or 'good lattice' multidimensional integration methods |35|-CT; and (b) 
integrations over plaques are accomplished by Monte Carlo methods. This results in the 
computational method that is simple and general, that lends itself to demonstration of 
numerical convergence, and that permits a systematic exploitation of coarse grained results. 

The a priori theoretical support for these models is limited; it has been mentioned 
above that they are physical and have a known simple connection to molecular theory. An 
additional point that is important for molecular justification of the continuum theory is 
that the interactions being treated, electrostatic interactions, are long-ranged on the scale 
of molecular sizes PT],H2[. Because the theoretical justification of these models is limited, 
empirical validation for the anticipated applications is important. Therefore, this paper 
presents a number of examples of dielectric model predictions of potentials of mean forces for 
rearrangements of molecular species and reactive complexes in water. It deserves emphasis 
that the issue of the validity of these models for predictions of changes in solvation free energy 
with changes in molecular geometry is distinct from those of accuracy for one geometry or 
of accuracy for changes in charge distribution at a fixed geometry. Most of the previous 
applications, electronic structure calculations and electron transfer theories, appear to rely 
on answers to the latter question. Preliminary results for several of the examples below have 
appeared in our previous work [p6| , p7[ . But with the exception of P3|Ji^ this is the only 
available work that obtains and studies potentials of mean forces, i. e. variations of solvation 
free energies with geometry, accurate on the thermal energy scale. Note also . Reference 
provided an important motivation for the present work. 
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2. Physical Model 

We consider a molecule in a solution environment and identify that molecule as the 
'solute' of interest. A solute volume is defined on the basis of its geometry. Partial charges 
describing the solute electric charge distribution are positioned within this volume. The 
Poisson equation of interest is 



V»£:(r)V$(r) 



47rp(r) 



(1) 



where p(r) is the density of electric charge associated with the solute molecule, the function 
e{r) gives the local value of the dielectric constant, and the solution $(r) is the electric 
potential. 



2.1. Integral equations for the model 

When contributions from an external surface vanish, the fundamental solutions of this 
Poisson equation are solutions of the integral equations: 



e(r)G(r,r') = G^''\r,r') + J V'G^^\r,r") 



e{r")G{r", v')dh\ 



(2a) 



G'(r,r') = G'(°)(r,r') - / \V'G^'^\r,r") 



'e{Y") - r 

An 



[V'G{r'\ r')] rfV", 



(2b) 



G(r,r')£(r') = G'(°)(r,r') + j G'(°)(r,r") 



y"e{v") 
47r£(r") 



V'G{T"y)e{j')dh". 



(2c) 



The derivation of these equations is discussed in an appendix. Here G(r, r') is the Green 
function for Eq. (|l|) — the electric potential at r due to a unit source at r' — and G'*^°^(r, r') 
is the Green function for the case where e{v) = 1. These equations apply both to the 
'unbounded case' and to periodic boundary conditions. For the latter case the potential is 
only required in a cell of volume V and the integral need only cover that volume. But in 
the latter case the electric potential due to a unit source "G*^'')(r, r')" is then the traditional 
Ewald potential. That result also is discussed in an appendix. It will be worthwhile to note 



that Eqs. (Q) are specializations of the more general relations that arise from considering 
how to calculate G'(r,r') for a specified e{v) with the help of a separately known G'^'^^(r, r') 
corresponding to £:*^°)(r). Those more general relations are: 

£(r)G(r,r') =£(°nr)G'(°)(r,r') 

+ / .(°)(r)V"G(°)(r,r") . (^)[V"ln (^-|i^)].(r")G(r", r')rfV", (3a) 



G(r, r') = G(°) (r, r') - J [v"G(°) (r, r")] • V'G{r", v')d'r\ (3b) 

V ^ ^ 

G(r,r')£(r') = G(°)(r,r')£(°nr') 

+ j^G^'\vy)e^'\v''){^)V'\n[0^^ (3c) 

Again, the derivation of these equations is discussed in an appendix. These equations will 
be of subsequent help in the design of numerical methods to treat systems for which the 
solute volume is not efficiently described as a finite collection of finite spheres, e.g., systems 
that include slabs and cylinders. 

2.2. Molecular Volume and Boundary 

In the customary applications of dielectric models a molecular volume is defined and 
the region outside this molecular volume is assigned a local value of the dielectric constant, 
e{r), equal to the measured dielectric constant of the solution. It is with the definition of 
the molecular volume that molecular scale information of the solvation structure, otherwise 
external to the model, is incorporated. It should be recognized that the defined molec- 
ular volume typically will change with the thermodynamic state of the system, with the 
structure of the solute, with solute-solvent interactions, and with correlations among sol- 
vent molecules. We comment later on the statistical thermodynamical underpinning of the 
model that clarifies the molecular origin of that information. The description of the molecu- 
lar surface should be simple and should permit disconnection with molecular rearrangement. 



Smoothness is desirable. The most popular previous choice has been the Connolly surface 
P6| defined as a surface of contact between a spherical probe and the van der Waals surface 
of the molecule. We have chosen instead to define the volume of a molecular complex in 
terms of a collection of spheres, each member s of that collection having an individually 
specified radius, R{s), and center, c(s). The volume enclosed by spheres in this collection 
is the 'molecular volume.' It is expected that the molecular volume will contain the source 
charges but that limitation could be straightforwardly relaxed. The molecular surface is 
defined as the boundary of the molecular volume. Generally, some of the surfaces of the 
individual spheres will be inside other spheres; some spherical surface will be buried and 
will not contribute to the molecular surface. 

Our considerations for making this choice of the molecular surface rather than the Con- 
nolly surface are that it is simpler and no less physical. The Connolly surface has the 
advantage of better smoothness for compact molecular conformations. But it has the dis- 
advantage that disconnection upon continuous conformational change either leaves sharp 
'cusps' or, if those cusps are arbitrarily shaved-off, the volume changes discontinuously. 
Since the statistical thermodynamic connection of the defined molecular surface to the sol- 
vation structure is somewhat indirect, involving solvent molecular sizes also, the differences 
between the present choice of molecular surface and the Connolly surface is not significant 
physically. 

This choice, and the Connolly choice also, yields a molecular surface that is generally 
not a Liapunov surface [0. That smoothness condition is an important feature of general 
abstract discussions of boundary element methods . The importance of the smoothness 
of the surface for general discussions can be exemplified by considering the case of the tradi- 
tional Dirichlet problem: solution of V'^ip{r) = within a closed volume with ip{r) specified 
on the boundary. If, for example, the boundary does not have a unique surface normal then 
discussions involving the normal derivative of ^p{r) on the boundary become vexed. The 
significance of this point for our applications can be discussed on the basis of Eq. ([2a|). The 



homogeneous contribution (7*^°^ (r, r') is continuous, smooth, and differentiable away from 
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sources and, on the assumption that no sources are located precisely on the surface, on the 
boundary of the molecular volume also. In order to isolate potential problems we can focus 
on the integral of Eq. pa|). We notice that the discontinuous function 6{r) appearing in 
the denominator can be cancelled exactly and that the factor G{r", r') is expected on the 
basis of Eq. (|^) to be continuous and finite away from sources and thus on the boundary in 
our problem. [In fact, we can typically arrange it to be non-negative. Then Eq. (^a]) can 
be a natural basis of Monte Carlo calculations.] If we agree not to require computation of 
the electric field precisely at points where no unique surface normal exists, then the only 
potential difficulty to consider is associated with the integrability of the surface dipole con- 
tribution to the electric potential at a point off the boundary at positions of the sources. 
But those contributions are clearly integrable and thus only a problem of practical difficulty. 

2.3. Comments 

Some points of note about Eqs. (0) are: first, Eq. ( pb| ) is simple to interpret physically. 
It says that the net electrostatic potential of a charge is composed of the potential due to 
polarization induced in the medium by electric fields there in addition to a direct, or bare 
interaction. Second, for the conventional applications 6{r) is a sharp step at the molecular 
surface and thus the integrals of Eqs. and (|2^) are surface integrals over the molecular 
surface. Although these equations are obviously twins, they offer different physical pictures 
of the problem. The integral of Eq. (|2a| ) describes a discontinuous function of r; that discon- 
tinuity is consistent with the obvious discontinuity of the left-side of Eq. (pa]). The integrand 
is a continuous function, however, for r not on the boundary. For r inside the molecular 
volume where e{r) = 1, the integrand of Eq. ( paf ) can be viewed as the potential due to a 
dipolar surface density G{r,r'). The integrand of Eq. (|2^) can be considered to be the po- 
tential due a surface charge density (47r£:(r"))^^ V"£:(r") • V"G{r", r')6{r'). For this case, the 
integrand varies discontinuously with normal displacement through the molecular surface 
because the normal electric field is discontinuous at a sharp step in e{r). If we evaluate that 
negative gradient normal to the molecular surface by assuming that it is half the sum of 
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the inside and outside values, we obtain from Eq. (p^) the equation that provides the basis 
for the most widely used of previous boundary element calculations p|,p|,P, p!0| , p^ , p^ , p!5| , p!8 



This is discussed in an appendix. Eq. can be written more physically for the electric 
potential of Eq. (0) as 

$(r) = (r) + J (r, r') (S^) • V'<l>(r')rfV' (4) 



V 



with the requirement that e{r) = 1 wherever p/(r) ^ 0. Here <I>*^°^ is the electric potential in 
the absense of the medium. 

In constructing and analyzing practical solutions is it helpful to recognize simple appli- 
cations of Gauss's law, for example to Eq. 



J V'G(°) (r, r') • V'e{r')d^r' = -47rA£r/(r) (5) 



V 



where ?7(r) is one or zero if the observation point (r) is inside or outside, respectively, the 
molecular volume. 

3. Examples 

All examples here used e = 77 A in the solvent. The quantities presented in these 
examples are the potentials of the mean forces obtained as 

W = U+{^) J p(r) ($(r) - $(0) {r))dh. (6) 

U is the static energy of the solute in the absence of the solvent. Typically, a constant 
contribution to W is adjusted so that this potential of mean force is zero at some asymptotic 
separation of interest. The distribution of free charge p(r) is a sum of partial charges so the 
integral of Eq. (|^) reduces to a sum over those charges. That is not an important limitation, 
clearly. Except where noted, we used the same charge distribution employed in obtaining 
the simulation results used as a comparison. 

All our results below were obtained on the basis of either Eqs. (Pi]), or (^), or both, 
utilizing methods described subsequently. Wherever possible we used the radii-parameters 
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of Ref . . Otherwise we adopted van der Waals radii arbitrarily and will note those values 



in the specification of the example calculations. We note that knowledge of the helpful 
work Ref. ||^ would probably have permitted better initial estimates of the required radii- 
parameters. However, we reiterate our expectation that variation of hydration free energies 
with geometry place a different demand on the dielectric model and upon radii-parameters 
than do total free energies or total free energy changes upon reaction. We expect that 
qualitative conclusions drawn below would not be changed by reasonable alterations of the 
radii-parameters used. 

Energy units on the following graphs are either k^T or kcal/mol. That choice was 
governed by the energy units chosen in reporting the original results that we are using for 
comparison. For all the results here k^T ^ 0.59kcal/mol and that smaller unit was typically 
adopted where more refined analyses were pursued. 

We note that practioners are not yet unanimous regarding the most appropriate boundary 
conditions to use for electrostatic interactions in simulation calculations of the sort we rely 
on here. Different boundary condition choices can lead to quite different results. We have 
exercised some judgement in choosing simulation results to use. In some cases, we have not 
used simulation results that are substantially the same as those we have used. In some other 
cases, however, we have excluded simulation results that display obvious dubious features, 
even though those features do not presently have a clear explanation. 

Our results below are based upon solutions of Eq. (P that are beheved to be accurate to 
about 0.3 k^T for the total solvation free energy for the worst cases. The principal difficulty 
we have encountered is the following: Numerical solution of Eq. (P requires replacement of 
a continuous problem with a discrete one. The results of different calculations can change 
practically discretely p3|-p5[1 unless special care is taken to avoid that and this is especially 



important for changes in geometry. For example, the position or number of lattice points 
can change discontinuously depending on the geometry and how the numerical problem is 
organized. Such changes can be small on the scale of the total solvation free energy, e.g. 
0.5% of the whole, but large and puzzling on the scale of the potentials of mean force. None 



of those features are evident in our results below. Thus, we believe our results presented 
here are the correct solution of the dielectric model for the problems considered. A major 
emphasis of the subsequent discussion of numerical methods is placed upon getting correct 
and correctly continuous variation of the free energy with geometry. 

3.1. Na+ • ■ • CI^ pair potential of mean force in water 



A variety of theoretical results [2(:,27,43,3Q-p3[ for the potential of mean force between 



a Na"*" ■ ■ ■ Cl^ ion pair in water are shown in Fig. |l]. We note that the XRISM ||5^ and the 
dielectric model calculations agree in properly describing the asymptotic 1/er behavior. All 
results show a contact minimum at about 2.6-2.8 A. The repulsive, short-distance side of 
this contact minimum is due to the effect of the solvent in pulling the ion pair apart. The 
direct ion core overlap repulsion begins to be significant only for distances less than 2.5 A. 

The two dielectric model calculations used different molecular surfaces. The earlier work 
of Rashin |Q employed a Connolly surface The later calculation used the van 



der Waals surface. The comparison of Fig. |l| shows that those differences in the molecular 
surfaces are important in the region of the free energy barrier. However, the close equality 
of those results near the contact minimum suggests that the differences in the molecular 
surfaces used are not important in that region of inter-ionic separations. 

On this basis, we restrict our conclusions here to the region of the contact minimum. We 
conclude that the dielectric models predict that contact-paired ions are too stable relative to 
the dissociation limit. The XRISM approximation is somewhat better in this regard but still 
errs in predicting too much stability for the contact ion pair. These errors are appreciable on 
the kT energy scale so that predictions of probability densities would be in error by factors 
of about 3. 

We note that including a modeled polarizability of the solvent and solute species does 
not seem to have a big effect on the predicted Na"*" ■ ■ ■ Cl~ pmf for distances less than about 

8A m. 
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3.2. CI ■ CI pair potential of mean force in water 

A variety of theoretical results 1 43. 5^, 52, 5^ of the potential of mean force between a 



C1+ ■ ■ ■ Cl^ ion pair in water are shown in Fig. ^ All results except 'molecular dynamics 
92' agree in properly describing the asymptotic 1/er behavior. Those exceptional molecular 
dynamics results were obtained at a non-zero salt concentration. 

The two dielectric model calculations used different molecular surfaces: Ref. |H3| adopted 



a Connolly surface [|^; the present adopted the van der Waals surface with the same radii- 
parameters. The dielectric model predictions are entirely repulsive. 

3.3. t-butyl+ ■ • ■ CI potential of mean force in water 

Fig. ^ shows theoretical results for approach of a Cl~ ion to a planar t-butyl"^ for ion- 
pairing configurations. The Monte Carlo and XRISM results are taken from Ref. For 



the dielectric model calculations, the radius for the tertiary carbon was taken as -R = l.SSA; 



for the united atom methyl groups those radii were R = 2.00A ^0|; for the chloride the 
Rashin-Honig [||] value R = 1.937A was used. Again, the dielectric model predicts too much 
stability of contact pairs relative to the dissociated fragments. The XRISM approximation 
is somewhat better but errs in the same direction by about 4kBT or about a factor of 55 in 
the associated radial distribution function. 



The results of Ref. P^ , Figure 5] cannot be conveniently placed on this graph. 



3.4. Na+ ■ ■ ■ dimethylphosphate ion pair potentials of mean force in water 

For the Na+ ■ ■ ■ dimethylphosphate" complex ion pair two sets of theoretical results 
are available |]56| , |57| . For both cases, the conformation of the dimethylphosphate anion 
was held rigidly in the gg structure of the fragment as described by electronic structure 
calculations |Q. In the first 'asymmetric' case the Na"*" ion was brought up along a P-0 



bond of an exposed oxygen atom; see Fig. ^. In the second 'symmetric' case the Na"*" ion was 
brought up along the bisector of the 0-P-O angle involving exposed oxygen atoms; see Fig. |^. 
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Beyond the sodium radius ||4^ of i? = 1.68A, the radii used in this case were van der Waals 
radii of Ref. phophorus, R = 2.10A; oxygen (ester), R = 1.65A; oxygen(phosphate), 



R = I.60A; united atom methyl, R = 2.00A. In the asymmetric case, the gaussian bath 
approximate theory agrees closely with the molecular dynamics results. The XRISM and 
dielectric model over-stabilize the asymmetric contact ion pairs. For the symmetric approach 
the molecular dynamics results are different from all the other theoretical results. Those 
theoretical results err by about 4kBT or about a factor of 55 in the first peak for the 
associated radial distribution function. 



We have not used the results of Ref. 



3.5. Chloride exchange in methyl chloride by symmetric SAr2 in water 

Fig. ^ shows available results for the model SAr2 reaction by a symmetric exchange of 
chloride in methylchloride along a linear reaction path [|6T| -|63[|. Beyond the chloride radius 



R = 1.937A [||], we used radii H: hydrogen, R = l.OOA; carbon, R = 1.85A. We used 



the partial charges of except for reaction coordinates in the range -lA< r < lA where 
used the revised model charges of Our dielectric model result appears to be in good 



agreement with the earlier work of Ref. |^ that was based upon AMI electronic structure 
model coupled to the dielectric medium. The agreement seen between dielectric models and 
the Monte Carlo results is good, better than that agreement between XRISM and Monte 
Carlo. 

3.6. Nucleophilic addition of hydroxide to formaldehyde in water 

Fig. shows available results for the attack of formaldehyde by hydroxide through a 
S7v2 mechanism ||6^j6^. We used radii hydrogen, R = l.OOA; oxygen (hydroxide). 



R = 1.65A; oxygen (carbonyl), R = I.6OA; carbon, R = 1.85A. The prediction of the 
dielectric model is qualitatively correct but with large quantitative errors. It has been 
shown that these quantitative errors can be empirically corrected by a reparameterization 
that pays no attention to physicality [^]. The XRISM results are quantitative in bettter 
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agreement with the Monte Carlo data but that XRISM approximation still incurs substantial 
quantitative errors, over-stabilizing the compact configuration relative to the dissociated 
species. 

3.7. Conformational equilibrium of N-methylacetamide in water 

All the examples here have been considered in order to test the theoretical model and 
not as a study of the particular system. Thus, the fidelity of the interaction model used 
for the molecular simulations to laboratory experiments has been less a concern than the 
correctness of the simulation results for the interaction model used. This point deserves 
particular emphasis for the N-methylacetamide example considered next . The conformations 
and interaction models that have been studied in this case appear not to be very realistic. 
However, the value of the simulation work as a benchmark result against which theoretical 
models can be compared remains. 

What has been done is to define two planes as follows: plane ^ 1 is defined as the 
plane of the H(N), CH3(N), and the C(carbonyl). Plane # 2 is the plane of the CH3(C), 
O, and the N. Then with the C(carbonyl) and the N atom fixed, all other bond lengths and 
angles fixed, and all force field parameters such as charges fixed, these planes are rotated 
relative to each other by an angle uj about the C(carbonyl)-N line. The trans conformation, 
which is the lowest energy conformation in the gas phase, corresponds to u = 180°. The cis 
conformation is = 0°. 

The infelicities of this study as a description of laboratory experiments are essentially 
two. The first was discussed in Ref. The description followed here using partial charges 
that are independent of conformation predicts that the trans isomer has a larger dipole 
moment than the cis isomer. Electronic structure calculations order the sizes these dipole 
moments oppositely [^,^. Because the modelled dipole moments do not display a correct 
dependence on conformation, the implied solvation contributions appear to be incorrect in 
preferentially stabilizing the trans conformation. When this feature of the molecular model 
is revised, the solvated free energies of the cis and trans conformers are substantially the 
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same 



The second difficulty has to do the with path chosen in these studies for interconversion 
of the two conformers of interest. A more detailed study would show that a natural path 
for interconversion of these isomers would be very different from the simple rotation of the 
planes defined above. Thus the results here should not be taken as a realistic description 
either of the net free energy difference between the two isomers or of the variation of free 
energy along a natural reaction path. 

However, a most important attribute of a test of theoretical model is that it is applied 
for precisely the same circumstances that characterize the experimental data. The available 
data are satisfactory for that purpose and such a comparison is given in Fig. ^ |6^. The 
radii-parameters used were Rh{n) = 1.25A, Racetyi-methyi = 1.955A, RN-methyi = 1-90A, 
Rc = 1.875A, R]^ = I.625A, Ro = 1.48A. The comparison shows that the dielectric model 
and all the other theories shown, except the XRISM approximation, give the right net 
solvation free energy. The apparent structural feature in the region 90° < u < 120° is 
missed by all the theories. 

4. Methods 



We specialize Eqs. (|2aD and (|2q ) or @ by taking the observation point r definitely 
just inside or outside that molecular surface defined above. In each case we then obtain a 
closed integral equation that can be solved to provide the information that determines the 
desired solution everywhere else. The integral in those boundary integral equations goes 
over the two dimensional surface of the molecular complex. The surface integral might be 
approximated by sampling the molecular surface using pseudo-random number series, quasi- 
random number series, or lattice rules These sampling points are equidistributed 
on the surface of the molecular volume. The development below eventually advocates using 
one of the latter two methods to define plaques on the molecular surface and then the use of 
Monte Carlo methods — with pseudo-random number series — to accomplish integrations 
over the plaques. For clarity of exposition, however, we ffist describe the direct 'sampling' of 
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the integrand. This samphng permits the evaluation of integrals over the molecular surface. 
In any of the cases anticipated, the integrand linearly involves the unknown quantities such 
as G{r,r'). Therefore, these methods of approximating the integral produce a finite set of 
linear equations that can be solved by standard methods. 

4.1. Direct Sampling of the Molecular Surface 

We begin by writing out the equations used with direct, or elemental, sampling of the 



molecular surface, treating Eq. (pa]) first. We introduce the notation 



/ G'(°)(r,r')p(r')rfV = <l>W(r), (7) 

J 

and similarly for $(r). We here focus on r infinitismally outside the molecular surface. We 
will denote the external (solution) dielectric constant by Eg and the molecular volume will 
be assigned a dielectric constant Em- An appropriate discretized version of Eq. (Pa|) is then 



£,$(r,) = $W(r,) + ^«;,,$(r,) (8) 

j 

where the coefficients are obtained as 

_ (r, - r,) • nj (e^ - eJ f A'KR{s^fri{v^) \ . . 

Here R{sj) is the radius of the sphere Sj whose surface is sampled by the jth sampling point 
and hj is the unit vector normal to the molecular surface at that sampled point. f]{rj) is 
an indicator function for the non-buried surface; it is zero if the jth sampling point is on 
buried surface and one if it is not. M{sj) is the number of sampling points used for sphere 
Sj. The left-most factor of £:(r) of Eq. (|2a| ) is Es because we are utilizing this equation for 
r just outside the molecular surface. The diagonal terms in this matrix of coefficients are 
simply obtained by noting that the potential just outside a uniform sheet of dipoles is 2ttx 
(the surface density of moment). Thus we have 

Wjj = 27T X X ri{rj) = (e, - Em)v{rj)/2. (10) 

47r 
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With these specifications, Eq. (Q) can be solved directly to yield the potential on the molec- 
ular surface. It is clear that sampling points on buried surface do not contribute to the 
physical resolution of the potential at the molecular surface beyond entering into evaluation 
of the matrix elements through the normalization numbers M{sj). As long as this nor- 
malization is properly maintained, sampling points on buried surfaces can be preemptively 
deleted from the list of sampling points. 

Note how this procedure works when applied to the simplest problem: the molecule is a 
single sphere with a charge at the center; Eg = e and = 1- We know that the potential 
is the same everywhere on the molecular surface. If we sample the surface only at only one 
point {M{sj) = 1) we find that the potential at the sampled point is 

1 1 2 

X 



R e-wii {e + l)R' 
The correct answer is 1/ eR. This example suggests the slight modification of Eq. (^) that 
arises from a slightly more refined estimate of the solid angle subtended by the molecular 
surface element at the jth sampling point: 

w,, = {e, - ej{l - M{s,r'/^)r]{r^)/2. (11) 

This refinement can be derived by assuming that the jth plaque has area 47ri?(sj)^/M(sj) 
and that area is bounded by a circle, i. e. each plaque is a circular cap with the spatial extent 
determined by the average. With this modification and only one sampled point {M{sj) = 1), 
the present method produces the exact answer for the potential at the surface when applied 
to this simplest problem. With this modification, the simple 'one-point' calculation satisfies 
the row sum rule T,kWjk = that follows from the application of Gauss's law Eq. (j^). 

This method works straightforwardly. The continuum model results of Ref. show 



what can be achieved with this primitive approach. We note that the difference between 
Eq. (p!0| ) and (0) was observed to be important in those calculations. One limitation is that 
to get accurate solutions the matrix Wij sometimes must be large, e.g. 10^x10^. For our 
present interests, 'accurate' has been gauged by the results for the solvation free energy. For 
the case = 1 this is 

16 



A/^ = Q I P(r) (<f (r) - <l>(°)(r)) d'r, (12) 

V 

should be correct on the thermal energy scale ksT and should display reasonable variations 
on this energy scale with variations of the solute geometry. Other properties of the solutions 
of these equations can be more sensitive. The electric fields in the interior of the solute 
are examples of such properties. We close this section by considering another of those 
more sensitive properties but still a property of physical interest, a molecular polarizability 
modelled on dielectric idealizations. 

The method above offers a simple way to model the polarizability of a solute molecule, 
by assigning a dielectric constant 7^ 1 to the molecular volume [^. This possibility 
is helpful for the physical fidelity of these models but also provides a tougher test of the 
numerical methods used. To establish a value of Em on the basis of a prescribed information 
about the net polarizability a we proceed as follows. Choose = 1 and augment $'-°)(r) to 
include a uniform external electric field: 

$(o)(r) ^ $(0)(r) -r (13) 

We then use the methods here to determine $(r) far from the solute molecule. More specif- 
ically, we calculate 

where the dipole moment p is computed from a specified origin. We note that the equation 
to be solved, consider Eq. ( pa|) for example, is linear. Therefore a response to a particular 
contribution $('^^(r) can be computed independently of other contributions. Thus, we can 
calculate without consideration of the charge distribution of the solute. It is a property 
only of the parameter and the shape of the molecular volume. It is independent of a 
coordinate system origin. With the specifications here we solve 

£„$^(ri) = -ri • + ^ Wij^r,{rj) (15) 
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as above, is a unit vector along axis rj. We then form the desired modeled net polarizability 
from 

a^r, = (1 - £m) e^, • ^ (r, - c{si)) ^^%{ri). (16) 

The distance of the plaque point from the center of its sphere c(sj) arises because the 
surface dipoles are oriented along the outward pointing normal to the surface. As an example, 
consider a sphere of radius R, sample the surface at only one point r and use Eq. ([TT|). The 
result is 

a,, = (..-l)i?^(^)(^). (17) 

The exact result is 

a,, = ('^^) R%,. (18) 

The trace of the approximate result is in error by a factor of 3 /{em + 2) that is natural. But 
the approximate result is crudely non-spherical. This property thus provides a tougher test 
of the numerical solution than does the solvation free energy because this property depends 
more sensitively on the uniformity of the sampled points. This property is therefore helpful 
in discriminating different solutions, e.g., with different levels of sampling. 

4.2. More Thorough Sampling for Integration over Plaques 

One way to go beyond the method described above to obtain more accurate computations 
is to evaluate the matrix elements more accurately. This step in achieving higher accuracy 
does not increase the dimension of the set of linear equations that must be solved; solving 
the linear equations is the computational limitation of the method above. We note further 
that the integrands are singular and, for typical molecular surfaces considered, discontinuous. 
These latter qualities are expected to make these integrals tough candidates for higher-order 
integration schemes. 

The method described here defines plaques, as suggested above, and then, in order to 
simplify the problem, assumes that one of the integrand factors can be treated as slowly 
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varying over each plaque. Such an assumption is satisfactory when the numbers of plaques 
are large and, consequently, the spatial extent of each plaque is small. Thus, we can be sure 
that the numerical solution converges to an accurate result as the number of plaques tends 
to infinity, just as with the method above. However, we expect the rate of convergence of 
these subsequent methods to be superior. 

We begin the development by defining plaques. The plaques will be defined as the Voronoi 
polyhedra of the plaque points on the non-buried surface of the each sphere. Typically the 
boundaries of the plaques will include the boundaries of intersection of spheres. The plaque 
points might be constructed with one of the 'sampling' methods mentioned above. For 
the present development it is not required that the method of construction of the plaque 
points be one that might be used for direct estimation of the integrals. It is natural though to 
include uniformity of distribution of the plaque points plus uniformity of size and shape of the 
plaques, and compactness, as desiderata. As specific examples, we have satisfactorily used 
both quasi-random number series as a sampling for the spherical surface, and the Lubotzky- 
Phillips-Sarnak spherical lattice to construct plaque points. We will denote the set 

of plaque points by {a}. 

In addition to the plaque points {a} we will require another larger set of points that 
sample the molecular surface uniformly and will be used to accomplish integrations over 
the plaques. We have typically used Monte Carlo methods with pseudo-random numbers 
to generate this set of points that will be denoted by {i}. When the plaque points are 
constructed to be reasonably uniformly distributed on the surface — even though only a 
finite set of points are used — we would typically include those plaque points in the set {i}. 



The present method can then be described by consideration of Eq. (paj). We assume that 
the G{r", r') appearing in the integral there is accurately a constant on each plaque. The 
equations for the off-diagonal coefficient Wap that follow from this assumption are 

The set of points i that are within the plaque j3 is then denoted by {i G j3}. M{sp) is the 
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number of points, sampling points plus the plaque point, on the sphere that supports plaque 
(3. Formulae for the diagonal element Waa that follow from the assumption that G(r", r') is 
constant over a plaque are 



The derivation of this result is depicted in Fig. ^ and discussed in an appendix. As indicated 
in detail by Eq. (|ll|), this applies to the case that the plaque point a is exposed. In Eq. (^0|), 
is the angle between the surface normals at point a, the plaque point, and the sampled 
point i. That formula arranges to use sample points outside the plaque a to calculate the 
solid angle subtended by that plaque at that plaque point. The purpose of that arrangement 
is to reduce the variance of that Monte Carlo estimation. Note that all sample points, buried 
or not, on the surface of the sphere that supports plaque a should be used in the estimation 
but whether any sample point resides within or without plaque a depends on resolution of 
buried surface. Eqs. ([T9|) , (pOf) , and (|2l|) are remarkably similar to Eqs. (|) and and 
remarkably simple. 

4.3. Refinement for higher resolution 

The left subscript index of Wap indicates a point on the molecular surface where the 
electrostatic potential is sought. It is straightforward to evaluate these matrix elements for 
points additional to the plaque points. It is straighforward to evalute the bare potential 
at additional points also. Evaluation of those quantities permits a simple calculation of a 
refined approximate $(r) once a preliminary or coarse solution, call it $(r), is obtained. We 
'iterate once' Eq. (||) with the preliminary solution used on the right: 




(20) 



and 




(21) 




(22) 
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In pushing to higher resolution, it is particularly helpful to note the requirement J2 = 
0. This follows from Gauss's law applied to Eq. (pa]) when r is outside the molecular volume. 



4.4. Coarsened equations 

A closed set of coarsened equations can be obtained from Eqs. (|22|) by the following 
argument: In the calculation of the matrix element Wap it was assumed that the potential 
was constant over the plaque (3. It is therefore conceptually consistent to consider $(r^) 
to be the average of the potential over the plaque rather than the value of the potential at 
one central point, the plaque point r^. With this choice of ^(r) as the plaque average, the 



plaque average of the finer resolution Eq. ( ]22|) then gives a coarse-grained equation that is 
closed in these plaque averaged quantities: 

£,$(r„) = $^°Va) + E^"/3^(r/3)- (23) 

/3 



Use of this equation prior to the higher resolution refinement Eq. (|2^) has the advantage 
that the coarse solution will then already satisfy the refinement in the plaque mean. The 
magnitude of the correction due to higher resolution refinement is then reduced. 

There is an alternative to this procedure that is typically advantageous, sometimes com- 
putationally simpler though requiring more memory, but also conceptually more indirect. 
We arrive at this alternative procedure by the following argument: we formulate a high res- 
olution set of discrete linear equations using Eqs. ([T9|), (pO]), and (^). We then appeal to a 



coarser tesselation of the molecular surface and require that all the points of the fine lattice 
that reside on the same coarse plaque have the same potential. This requirement produces 
an overdetermined set of linear equations that can be analyzed with the traditional applica- 
tion of singular value decomposition |]TD[ to produce the solution that minimizes the mean 
square residual of fine resolution equations. This approach attempts to utilize the fine-scale 
spatial variation of the $*^°)(ra) in constructing a coarse solution. 
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4.5. Surface Charges 

An alternative implementation of Eq. (Hf) holds an advantage for applications of tradi- 



tional electronic structure packages ||7T|. That advantage derives from the fact that the 
induced electrostatic fields are represented with surface distributions of charge. Traditional 
electronic structure packages can use collections of charge monopoles simply. If the surface 
charge distributions are represented by a collection of monopoles then the utilization of tra- 
ditional electronic structure packages is logistically simplified. In this section we describe 
how the ideas above are adapted to utilization of Eq. in those applications. 

As with the development above we first specialize Eq. (^) to extract a boundary integral 
equation. This is discussed in further details in the appendix. We evaluate formally the 
gradient of that equation and then consider how to interpret the results for e{r) a sharp step 
at the molecule surface. We interpret integrand factors that become discontinuous in that 
circumstance as the average of the values inside and outside that discontinuity. We then 
take the observation point r just inside the molecular surface and notice that we obtain the 
simple closed equation 

n • D(r) = n • Eo(r) + / n • VGo(r, r') (^^^^) n • D(r')rfV (24) 

s 

for the Maxwell displacement normal to the molecular surface at the point r. E^{r) is 
the component of the bare electric field — due to p — normal to the molecular surface at 
the point r. Eq. (^) involves a surface integral only because it is already specialized to 
treat a sharp surface. As noted above, this equation provides the basis for the most widely 
used previous boundary element approaches for these problems [jI|J^J^,0,|l^,|l^,|^,0] . The 
discretized version of Eq. (plD, comparable to Eq. (H), is 

D„(r,) = E°(r,) + Y.^^JDr,{r,), (25) 

j 

where the off-diagonal coefficients are given by 

_ h, • (r, - r,) - ( A7TRis,fr]ir,) \ 



|rj — Fjl'^ Aire 
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The diagonal contributions can be obtained by noting that when the plaques are small, their 
curvature can be neglected and that, since the plaque point is immediately inside the plaque, 
edge effects may be neglected in a first approximation. The integration required in Eq. (|2^ ) is 
then minus the electric field due to a planar sheet with charge density {sg — Sm) Dnij) /Anes- 
This electric field just inside has magnitude 27rx (that surface charge density) and is anti- 
parallel to the outward directed surface normal. Thus 

W.jj = {Es - £rrMrj)/'2£s- (27) 

The refinement corresponding to Eq. (|lT]) is here 

wjj = {s, - eJil - M{s,)-^'^)r]{r^)/2es. (28) 

After Dn is obtained on the molecular surface, we review the antecedents of Eq. (p^ to see 
how to obtain the thermodynamic quantity A/i of Eq. (p!2D. When all the free charges are 
in regions where £:(r) = e^, the induced potential should be calculated as 




in close similarity to Eq. (^. Again, the refinement of Eq. ( pH]) produces the exact answer 
for the simplest problem, the spherical ion with charge at the center mentioned above, when 
the surface is sampled at one point only. 

4.6. Coarse and Finer Solutions with Surface Charges 

The arguments above that identified a higher resolution approximate and the 'self- 
consistent' coarse equations carry through for the present surface charge of the adjoint 
prespective also. 

5. Molecular theoretical perspective 

Here we note some points of prespective gained by an appreciation of how the dielectric 
models for solvation free energies connect to molecular theory. Some of these points have 
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been noted before p6|- p9| , pl| -|3^ ■ But they seem not be widely appreciated. Further, new 



points here clearly indicate how to improve simply the predictions of potentials of mean 
force while utilizing the methods detailed above. 

The molecular theory corresponding to these dielectric models is the second-order cumu- 
lant approximation 

A/.^A/.o + <^E^(j)) -f ((|E^a)-(E^(k)) j ) ■ (30) 

Here is the electrostatic interaction potential energy coupling the solute to solvent 
molecule j. The brackets (■ ■ ■)q indicate the thermal average in the absence of those elec- 
trostatic couplings and A/xq is the excess chemical potential of the solute at infinite dilution 
again in the absence of electrostatic interactions. When the solute is an infinitely dilute 
second component this molecular approximation is perturbation theory through second or- 
der in the electrostatic interactions. When the molecule is not literally an infinitely dilute 
second component the 'infinite dilution' restriction means that one molecule is distinguished 
for the purposes of calculation. This is still a natural theory but the medium now contains 
a non-zero concentration of molecules mechanically identical to the 'solute.' The medium 
properties non-perturbatively reflect that fact. From the perspective of this molecular theory, 
the dielectric solvation model application neglects the zeroth and first order terms, makes an 



estimate of the second-order term, and neglects all succeeding contributions. Refs. p8| , |29|j31 
can be consulted for alternative views. 

In view of the generally satisfactory predictions of the dielectric models in the examples 
above, it is reasonable to conclude that the approximate molecular theory of Eq. (|30|) would 



provide a valid description of those systems also but with further molecular detail p3|j34 



and without the empiricism encoded in the radii-parameters of the dielectric model. That 
empiricism does in fact limit the ability to use the model for prediction. In the examples 
of the inter-ionic potentials of mean force above, it is clear that the dielectric model is not 
satisfactorily accurate. But even then, it is a reasonable view that the failures of the model 
are due principally to failure in suitably describing the molecular volume rather than an 

24 



essential failure of the perturbative formula Eq. (|30|). 

A direct implementation of Eq. (|30| ) would require simulation calculation. Such calcula- 
tions are likely only to be somewhat simpler than the direct simulation of the potential of 
mean force sought. Thus, to extend the theoretical insight that Eq. (^) is physically valid, 
other tricks must be found. One additional idea is to attempt to calculate the first-order 
and second-order terms in the molecular theory on the basis of the dielectric models. Direct 
use of that idea is likely to be unsuccessful because those terms in the molecular theory 
require analysis of the solvation of a discharged solute and the dielectric model is unlikely 
to be accurate for discharged solutes. However, a physically equivalent but operationally 



distinct expression of that idea is 1 72,73 



A/.^ A/io + (E^(j)) ■ (31) 

\ i ' 1/2 

With this formula, the average can be evaluated by considering one solvent molecule, a probe 
molecule, in addition to the solute discharged to the extent 1/2. This molecular grouping 
can be studied on the basis of the dielectric model with the expectation that the dielectric 
model should give a physically valid description. The dielectric model can provide a pair 
potential of mean force for the additional solvent molecule relative to the half-discharged 
solute. The mean electrostatic potential that the water molecule expresses on the solute is 
then the quantity of interest. To obtain that quantity requires averaging over the Boltzmann 
factor of the dielectric model potential of mean force for the grouping. 

The justification of this approach is that the formula Eq. (|3TD would give precisely the 
same answer as Eq. (|30D if that latter result were precisely correct. However, if we are looking 
for one such formula to exploit with the dielectric model, Eq (^) is preferred because: (a) it 
uses the dielectric model in a region where it might be expected to be physically valid; and 
(b) Eq (^) requires considering only one water molecule in addition to the solute rather 
than two as would the last term of Eq (|30D. 

More formally, we can argue as follows. The desired excess chemical potential can be 
expressed as 
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A/. = A/io + I'lr^vii)/ d^- (32) 

C, G [0, 1] is a parameter that scales all the partial charges of the solute. The mean value 
theorem applied to the integral yields 

A/i = A/io + (E^(j)) , (33) 

where < ^* < 1. The success of the dielectric model suggests that we use the estimate 
^* ^ 1/2 and evaluate the average by exploiting the dielectric model itself. This more formal 
argument obviously further suggests that the precise value used for though near 1/2, can 
be helpfully empiricized as sufficient experience is acquired with these methods. 

In closing this section it is helpful to note some often mentioned alternatives for the 
molecular theory corresponding to this dielectric model. One alternative might be based 
upon the suggestion that a spatially local dielectric constant should be used. If a first 
solvation shell is structurally fixed and if a local dielectric constant is given to describe 
the polarization of that shell, then this is similar in spirit to the present model. However, 
more generally it must be admitted that a local dielectric constant is a ambiguous concept 
unless local means 'macroscopically local.' The macroscopic solution dielectric constant is 
unambiguous and the net polarization of the interior of a rigid molecular structure also is a 
simple physical concept. This level of modeling is therefore reasonable but we have refrained 
from extending the present model to elaborate those ideas in more difficult settings. 

A second alternative is based upon the observation that on a molecular scale the dielectric 
response should be non-local. This non-local response would be characterized by a wave- 
vector dependent dielectric constant, e{k). This is a much less ambiguous approach than 
the 'local dielectric constant' idea mentioned above. However, it does not go as far as 
the second-order perturbation theory advocated here. The 'non-local' response approach 
would not typically include the linear perturbative term. The non-local response approach 
would typically compromise the description of short-ranged 'packing' interactions as they 
are expressed in the perturbative terms. Moreover, the i{k) is not necessarily simple |]7^ , |75| . 
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6. Conclusions 

This dielectric model |jl|J^ gives interesting and helpful results for the variation in sol- 
vation free energy with solute geometry. However, it typically significantly over-stabilizes 
classic attractive ion-pairing configurations. 

In several cases, e.g. Fig. ^ and Fig. ||, the predictions of this dielectric model are quite 
similar to those of more ambitious theories. This comparison clearly teaches us about the 
performance of those more complicated theories for those cases: the phenomena described 
are unsuspectedly simple and those more complicated theories are at least qualitatively valid 
for those simple circumstances. 

It is numerically feasible to solve this dielectric model to thermal accuracies for the excess 
chemical potential pB|^5|. The integral equations Eqs. (0) and (^) that are the basis for the 
numerical boundary integral methods studied here have not been given previously in this 
context. The algorithmic discussion of Eq. (|^) has laid the basis for a multigrid method of 
solution. Such methods are likely to be necessary for similarly accurate numerical solution 
of these models for much larger solutes. Eq. @ will be particularly helpful for macroscopic 
surfaces such as solution interfaces and membranes. Eqs. (^8]), (|63D, and ( 16^) show how 
these methods can be transferred directly to periodic boundary conditions. That will be 
important for two reasons. First, the best controlled quantitative results on hydration and 
electrostatic interactions utilize periodic boundary conditions on electrostatic interactions. 
Second, simulation of the hydration of macromolecules often involves periodic boundary 
conditions, only a few layers of hydrating water, and almost no direct investigation of system 
size effects; treatment of hydration on the basis of the dielectric model with and without 
periodic boundary conditions would be a natural first step towards getting quantitative 
thermodynamic limiting results for such solutions. 

Finally, we draw attention to the connection discussed here of this dielectric model to 
molecular theory. As a proposal for going beyond current applications of dielectric models, 
the probe solvent molecule Eq. (pID is new. By reintroducing molecular solvation structure. 
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it is likely to correct the largest part of the errors of the dielectric model, i. e., the quantitative 
inaccuracy for either molecular surface used of the contact minimum in Fig. |1] and of the 
barrier to escape from that contact minimum. As a guide in analyzing molecular statistical 
thermodynamic results, this relation has already been useful. At a more practical level, 
this approach should reduce the sensitivity of the model results on the radii-parameters. 
It should be cautioned however, that the dielectric model is clearly approximate and the 
quadratic theory also is approximate [^|4l- Typically, these results can be expected to be 



superceded for applications where molecular accuracy is important. 
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Appendix A: Integral formulations 

We wish to write the equations to be solved as integral equations because those integral 
equations serve as the basis of the numerical methods here and because this is the most 
natural way to make contact with molecular theory. We consider the equation for the Green 
function, ^(r, r'): 

V«£(r)VG(r,r') = -47r(5(r-r'). (34) 

G{r, r') is the potential at r due to a unit charge at r'. We can rearrange this into a 
perturbation theory standardly. Consider the reference problem 

eoV • VG(°)(r, r') = -47r5(r - r'). (35) 
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G^^\r, r') is the potential at r due to a unit charge at r' for a uniform dielectric with dielectric 
constant Eq. Then the perturbation theory is 



^(r, r') = G^o^r, r') + J G^''\r, r")V" ■ 

V 



sir" 



47r 



Y'G{r", r') 



'3/1 



(36) 



When Sq — 1 we interpret this equation by noting that 



V". 



47r 



(37) 



is then the induced charge at r" due to a unit charge at r'. The total potential ^(r, r') at r 
is the sum of the bare potential G^^\v, r') of the unit source charge plus the bare potential 
of all the induced charges. 

An integration by parts casts this into the form: 



^(r, r') = G^°)(r, r') - j [VG^^^r, r")] • ( ^^''4^ [V"G(r", r')] d?r" . 

V ^ ^ 



(38) 



The surface contributions vanish because we are here considering a charge distribution 
of finite extent and the potential is required to vanish at infinity. This form makes evident 
the well-known symmetry of G(r, r'). Now when Eq = 1 the interpretation is that 

(e{v")-e,' 



47r 



[-V"G(r",rO] 



(39) 



is the induced polarization — 47rP = {s — 1)E. Since V"G^^\y, r") • u is the potential at r 
due to a unit dipole u at r", the total potential G(r, r') at r is the sum of the bare potential 
(^{o)(r^ r) of the unit source charge plus the bare potential of all the induced polarization. 
The form 



^(r, rO = ^("^(r, r') + j |v" • (^^^^^^) ^^^(r, r"; 



G{r",r')dh" 



(40) 



is obtained by another integration by parts. By elaborating the gradient operations we see 
that 



V". 



e r" 



£0 



47r 



V"G'('^)(r,r") 



V" 



47r 



V"G'('^)(r,r")- 



S{r-r"). 
(41) 
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Our basic equation then becomes 



V 



• 


"V"£(r")' 




_47r£(r")_ 



[e{r")G{r", r')] d^r". (42) 



eoG^^\v, r') is the bare Green function for a uniform dielectric with dielectric constant one. 
Thus we write 



e{r)G{r,r') = G'~^\r,r') + J V'G'^^\r,r") 



e{Y")G{v\ v')dh\ 



(43) 



with the understanding that here G'^'^)(r,r') is the vacuum Green function. This integral 
form of Poisson's equation directly gives the correct answer for a uniform dielectric. It also 
admits a simple interpretation. [V"£:(r")/47r] G(r", r') can be viewed as a surface density 
of dipoles. The remaining integrand factor W'G^^^ (r, r") is then the potential due to those 
surface dipoles. This formula has the following consistency: the indirect term on the right- 
side, that associated with the surface dipole distribution, is discontinuous across the surface. 
The value of such a discontinuity is 47rx (the local surface density of dipoles). Here that is 
47r[(£: — l)/(47r)]G'(r, r') for r located on the surface. The left-side of this equation also is 
discontinuous and the amount of that discontinuity is {e — l)G(r, r'). 
If we manipulate the Eq. (^) this way we obtain 

V" . Mr") - So) V'G{t", t')] = [V"e(r")] • [V"G(r", r')] 

+ ( '^"""h.'' ] {-47r5(r" - r') - [V"e{r")] . [V"G(r", r')]} . (44) 



e{r") 



Thus 



G(r,r')e(r') = G^°\r,r') + J G'-^\r,r") 



47r£(r") 



V'G{r", r')eir')d\". 



(45) 



Here again we understand that G^'^\r,r') is the vacuum Green function. This equa- 
tion can be viewed as the basis of the previous boundary element numerical approaches 
||l|j3|J^,p!0|, |T^ , |T^ , p!5| , p!^ . If the sources are located at positions r' where the dielectric function 
e{r') is one then this equation states that the full potential can be composed as the bare 
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potential of those sources plus the bare potential due to charges distributed throughout 
regions of inhomogeneity of e{r). 

To show the explicit correspondence with previous workers, we assume that all charges 
are located in regions for which £(r) = 1, right-sum over those charges, operate on Eq. (12) 
with — fi • V for r just inside the molecular surface to obtain 

'e-V 



n • D(r) = n • Eo(r) + / n • VG^°\r, r') (^^) n • D(r')dV'. 



(46) 



The right-most quantity in the integrand — the surface density of charge — is worked out 



as: 



?(r') 



e-1 



An J \£ + lJ '2 
^^')n«D(r') 



^ ^ 'J)(n«D(r')/£ + n«D(r'))) 



(47) 



These integral formulations can be generalized to cases in which the zeroth order contri- 
bution is a solution for a similar problem with nonconstant £{r). Such generalizations will 
be helpful in treating molecules near surfaces. The required developments here follow the 
above work closely. We reconsider the reference problem and select 



V • £o(r)VG(°Hr,r') = -47r5(r - r'). 



(48) 



G^^\r,r') is the potential at r due to a unit charge at r' for the reference problem with 
dielectric constant eo{r). Then the perturbation theory is 

'Ae{r") 

tr^ ' l^r, r ; V • i 

V 



Gir,r') = G'(°)(r,r') + J ^(^^(r, r")V" 



47r 



^^(r^r') 



(49) 



where we use the notation A£:(r) = e{r) — eo{r). We then proceed in the obvious way: 



^(r, r') = G^o^r, r') + / G^o^r, r")V" • \( Mp\ £{r")V'G(r", r') 

i [\4:TT£{r")J 



(50) 



Working out the V's gives 



G(r,rO = G(°)(r,r')£o(r')/£(rO +/ G^'\r,r")V' (|^) . e(r")V"G(r", rO^V. (51) 
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But 

Then finally 

G(r,r')5(rO=G(°)(r,r')eo(r') 

+ / GW(r,r")eo(r")V"ln (^].V'G{r"y)e{v')£r". (53) 



V 



This is the desired generalization of Eq. (|45| ) . 

For the generalization of Eq. (^31), we proceed similarly but start from the adjoint form: 



G(r,r')=G(°)(r,r')+/ |v" . [(^^) V"G(°)(r, r")] } G(r", r' 



)ciV" (54) 



We then find 



£(r)G(r,rO=eo(r)G(°)(r,r') 

+ J eo(r)V"G(°)(r,r") . (^)[V"ln (^^^]e(r")G(r", r')t/V'. (55) 

Appendix B: Boundary integrals, periodic boundary conditions, and dielec- 
tric models 

The methods above can be used also to solve these dielectric models in periodic boundary 
conditions. This may be important because almost all simulations of macroscopic matter 
utilize periodic boundary conditions; and many of those use the solution of Poisson's equation 
in periodic boundary conditions, the Ewald potential. 

We start by noting that the Ewald potential could be constructed with a boundary 
integral method. This can be demonstrated in a classic fashion from Green's second formula 

El 

J rfV{/(r)V2(7(r)-(7(r)VV(r)} = / rfV n • {/(r) V(7(r) - (?(r) V/(r)} (56) 

V dV 

in the present notation. The integral on the right is over the surface of a unit cell of a Bravais 
lattice, e.g., the surface of a cube, n is the outward pointing normal. We exploit this in the 

32 



traditional fashion with the traditional choice /(r) = l/|r — r'| and g{r) = '^Ewaidi'^)- Then 
V^/(r) = — 47r5(r — r') and V'^gir) = —An (5(r — r') — V^^). This gives the equation 

^E^aidir) = -- T7T——ii + drn.l ^En^aidir )V — — . (57) 



r 

V ' ' dV V 1 1 11^ 

This notation assumes, and our argument depends on this construction, that the unit cell is 

centered on the source. This will not limit the generality of the result because we will use 

that device to construct the general periodic solution of interest. With this choice of origin, 

the Ewald electric field normal to the surface must be zero. Then 

1 r d'^r' f 1 

VEwaid{r) = / — - / d^r' h»(fEwaidir')V— (58) 

r J V \r — r\ J r — r 

V ' ' dv ' ' 

Thus, the Ewald potential can be represented as the Coulomb potential due to a charge 
distribution within the cell plus the Coulomb potential due to a dipolar surface distribution. 
This equation might be used to determine the Ewald potential by discretizing the surface, 
taking the position r on the surface, and using Eq. (|58|) to determine ^EwaidiT") on the 
surface. The periodic solution could then be constructed by periodic replication. However, 
this equation is properly deficient in one respect: it is ambiguous with respect to the addition 
of a spatial constant. Thus if ^(r) is a solution, so is ^(r) + C. It is traditional |]78[ and 
natural |7S| to chose this constant so that 

/ ^Ev.aid{r)d^r = 0. (59) 

V 

With this backdrop, we reconsider the dielectric problem. A relevant Green's formula is 



now 



J d\ {/(r)V . W)Vg{v)] - g{v)V • [^(r) V/(r)]} 

V 

= Jd\h. {f{r)e{r)Vg{r) - g{r)e{r)V f{r)} (60) 



We choose 



/(r) = VEn^aidir - r'), V'f{r) = -Atc (5(r - r') - 

g{r) = <^(r, r"), V • e{r)Vg{r) = -4n (6{r - r") - V-') . (61) 
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The final point to make about this application of Green's formulae is that in this case the 
surface contributions will vanish: 

J d\h* {fEwaidir - r')e{r)Vv{r, r") - ^(r, Y")e{Y)V ^En>au{^ " 0} = 0- (62) 

dV 

The reason for this is that all functions in the integrand are triply-periodic. Thus, the 
gradient at a surface point can always be matched against another with equal weight and 
opposite projection on the surface normal. The relations V»e(r)V/(r) = (V£:(r))«(V/(r)) + 
5(r)VV(r) and V5£;«,aZd(r - r') = ipEwaid{r' - r) then directly lead to 

eivMr, r') = ^.^^...(r - r') + J d'r"V' ^Eu^aui^ - r") • V'e{j")^{j'\ r') 

V 

+ V-^ J [e{v)if{v, r') - ipE.naid{v - r')] . (63) 

V 

This should be compared to Eq. (|43|) . Now if we use the conventional Ewald potential, 
then we have Eq. (|5^). Further, if we use Eq. (^) without the penultimate contribution, 
then e{v')ip{v') is a linear combination of just those conventional potentials. In that case, 
the neglected last term will be zero. Our conclusion is that the form of the equation need 
not change in going from the unbounded to the Ewald case provided we substitute the 
conventional Ewald potential for Coulomb's law. 

A similar calculation starting from Eq. (p6D yields 

¥.(r, r') = ^E^au{r - v')/e{v') + j d^r'^E^auir - r") {^^^ • VV(r", r') 

+ I d\ [(^(r, r") - ^E^aidiv - v")/e{v)] . (64) 

V 

Again, if we use the conventional Ewald potential, the last term on the right vanishes. 
Appendix C: Solid angle from sampled surface 

Here we consider the formula Eq. (^). Our goal is to show how to use points sampled 
uniformly on the surface of a sphere to estimate the solid angle subtended by a plaque at 
the plaque point just outside the sphere. 
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Fig. ^ shows the geometry considered here. We take the z-axis through the plaque point 
with origin at the center of the sphere. The comphcation is that the sohd angle required 
is referenced to the origin at the plaque point. The solid angle subtended just inside the 
plaque point is given by 

n.ns.,e = 2n + jHiu)du (65) 

where du is denotes the integration over the surface of a unit sphere centered at the plaque 
point and H{uj) is the indicator function for the angles u for which the unit vector u points 
within the plaque. This integration can be implemented as 

/•27r /-O 

^inside = 2n + dip J ^d {cos '&)H{(^,'d). (66) 

In order to use sampling points uniformly distributed on the surface of the sphere we need 
to reference the integration to the center of the sphere rather than the plaque point. In that 
change of angular coordinates the azimuthal angle if is unchanged but the polar angle 
transforms as 



d (cos ^) = d{cos^s)/ p^/yi - cos^. 
The polar angle from the center of the sphere is i)s- We then write 

r2TT rO 



(67) 



pZlT pU 

^inside = 2^ + J I ^"^^ (C0S^9) 



C J\ H(^, m (cos m 

J^"" d^j\d {cos ^ ' 

This is in a form of an average: 

^inside = 27r + ^ ((1 - cos^9,)-'/')^ (69) 

where we have used the transformation Eq. (0). The subscript indicates that the sampling 
points should be uniformly distributed on the surface of the sphere. 

It can be checked that this formula produces the correct answer when the plaque covers 
the entire sphere. When the plaque is a circular cap centered at the plaque point this formula 
produces 
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^inside = 2n(^l + v^(l-cosi?*)/2) (70) 

where is the polar angle at the perimeter of the circular cap. The use of the engineering 
estimate 



'"xfil (71) 



(1 - cosr) /2 = 27r (1 - cos^) /An 

^^N' V47r 

underlies the approximation Eq. (|TT|). 

The integrand of Eq. (|69D is integrable but singular. This difficulty can be relieved by 
considering the complementary indicator function 1 — H{(f, ■(9) and the complementary solid 
angle: 

f27r rO 



'n-H{^,^)){l-cosAr'^')^. (72) 



23/2 

But this complement is just the outside solid angle desired, floutside = An — flmside- Imple- 
menting this formula as a sampling estimate that uses the points uniformly distributed on 
the spherical area outside the plaque gives Eq. (EOl). 
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FIGURES 



Figure 1. Theoretical results for tlie Na+ ■ ■ -Cl^ potential of mean force in water under 
physiological conditions. The two molecular dynamics results differ principally in the fact 
that the upper molecular dynamics curve was obtained at non-zero concentration of salt; 
see Ref. [Q. The curve labeled 'molecular dynamics 91' is the result of Ref. [^]. A simple 
translation nearly suffices to superpose these two molecular dynamics results. The XRISM 



curve is redrawn from Ref. ||53|. 'Dielectric 89' is redrawn from Ref. and 'dielectric 



94' is redrawn from Refs. p6| , p7| . The two dielectric calculations differ in that the earlier 
one used the Connolly molecular surface and the later one used the van der Waals surface. 
Those results are accurately the same in the region of the contact minimum but the different 
surfaces imply appreciably different results in the region of the free energy barrier. 



Figure 2. Theoretical results for the Cl~ ■ • -Cl^ potential of mean force in water under 
physiological conditions. The two molecular dynamics results differ principally in the fact 
that the lower molecular dynamics curve was obtained at non-zero concentration of salt 
(Ref. [^). The curve labeled 'molecular dynamics 91' is the result of Ref. [^. A simple 



translation nearly suffices to superpose these two molecular dynamics results. The XRISM 



curve is redrawn from Ref. |53|. 'Dielectric 89' is redrawn from Ref. and 'dielectric 
95' is new here. The two dielectric calculations differ in the molecular surface assumed as 
discussed in the text and the caption of Fig. |l]. 

Figure 3. Theoretical potentials of mean force for approach of a Cl^ ion to a planar 
t-butyl+ ion. See Refs. MEB- 
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Figure 4. Theoretical potentials of mean force for approach of Na"*" along the P-0 
bond of an exposed oxygen atom in methylphosphate" — the asymmetric approach. See 
Refs. |5ff8[. 



Figure 5. Theoretical potentials of mean force for approach of Na"*" along the bisector 
of the 0-P-O bond angle of exposed oxygen atoms in methylphosphate" — the symmetric 



approach. See Refs. ||56|-^ 



Figure 6. Theoretical potentials of mean force for exchange of CI in methylchloride by 



a symmetric S Ar2 process along a linear reaction path. See Refs. p7| , ^|6T| -|63[| . The reaction 
coordinate in the difference between the two carbon-chloride distances. 



Figure 7. Theoretical potentials of mean force for nucleophilic addition of hydroxide to 
formaldehyde. See Refs. P7| , |5^|53[ . The reaction coordinate plotted is the carbon-oxygen 
distance. 



Figure 8. Solvation contribution to the potential of mean forces for model interconversion 



of cis [left, LU = 0°] to trans [right, lu = 180°]. See text and Refs. [pq-pS 



Figure 9. Geometry for calculation of the solid angle subtended at a plaque point. 
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